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ABSTRACT 

Aims. In light of the recent detection of direct evidence for the formation of Kelvin-Helmholtz instabilities in the Orion nebula, we 
expand upon previous modelling efforts by numerically simulating the shear-flow driven gas and dust dynamics in locations where 
the H// region and the molecular cloud interact. We aim to directly confront the simulation results with the infrared observations. 
Methods. To numerically model the onset and full nonlinear development of the Kelvin-Helmholtz instability we take the setup 
proposed to interpret the observations, and adjust it to a full 3D hydrodynamical simulation that includes the dynamics of gas as well 
as dust. A dust grain distribution with sizes between 5-250 nm is used, exploiting the gas-hdust module of the MPI-AMRVAC code, 
in which the dust species are represented by several pressureless dust fluids. The evolution of the model is followed well into the 
nonlinear phase. The output of these simulations is then used as input for the SKIRT dust radiative transfer code to obtain infrared 
images at several stages of the evolution, which can be compared to the observations. 

Results. We confirm that a 3D Kelvin-Helmholtz instability is able to develop in the proposed setup, and that the formation of 
the instability is not inhibited by the addition of dust. Kelvin-Helmholtz billows form at the end of the linear phase, and synthetic 
observations of the billows show striking similarities to the infrared observations. It is pointed out that the high density dust regions 
preferentially collect on the flanks of the billows. To get agreement with the observed Kelvin-Helmholtz ripples, the assumed geometry 
between the background radiation, the billows and the observer is seen to be of critical importance. 
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1. Introduction 


Sometimes a little push is all that is needed to make a seemingly 
stable fluid evolve into a turbulent state. Typically this transition 
is caused by a fluid instability, and many of these mechanisms 
have been studied extensively in the past decades (see e.g. jChan^ 
drasekhar] ( |1961| )). The Kelvin-Helmholtz instability (KHI) is a 
notable example of this as it plays an important role in a wide 
range of different fluid applications such as for example oceanic 
circul ation ([van Haren & Gostiaux||2010| ), winds on planet sur¬ 
faces ( [Chapman & Browning 1997), the flanks of expanding 
coronal mass ejections ( jFoullon et al. 201 Ij), mag netic reconnec¬ 
tion in the solar corona ( Lapenta & Knoll|2003|), intera ction be¬ 
tween comet tails and the solar wind ( Ershkovich|1980j), mixing 
of solar wind material into Earth’s magnetosphere ( jHasegawaj 
et al.|2004] ), astrophysical jets ( jBaty & Keppens|2006| ) and many 
others. While the KHI is a hydrodynamical instability, magnetic 
fields can alter its dynamics and cause stabilisation or further 
destabilise the setup. As the previous range of examples demon¬ 
strates, many of the relevant astrophysical fluids in the KHI is 
of importance display magnetic effects. In molecular clouds, the 
KHI has been linked to the formation of filamentary structures 
( Hendrix & Keppens|2014| ), as well as to turbulence formation. 
While the source of turbulence, observed in molecular clouds 
through the detection of non-thermal line-widths around 1x10^ 

- 2 X 10^ cm s“\ is still debated, it has been linked at least par¬ 
tially to the KHI allowing to transfer energy to smaller scale 
structures ( jElmegreen & Scalo]|2004[ [Berne et aL]|2011t [Berne 


[& Matsumoto|2Q12] ). While the occurrence of the KHI in space 
is clearly established, direct evidence of ongo ing instabilities ar e 
harder to obtain. At a distance of 412 pc ( [Reid et ^[2009[ ), 
the Orion nebula is the closest Hjj region. Its association with 
young massive stars and its apparent brightness make it an in¬ 
tensively investigated region over a large range of frequencies 
( [Q’dell[[T001[ ). As such, it is an ideal laboratory for investiga¬ 
tion^ smaller scale structure development. Recently [Berne et aL] 
( [2010[ ) discussed mid-infrared observations of ripple-like struc¬ 
tures on the edge of the Orion nebula’s Hjj region and the sur¬ 
rounding giant molecular clouds. The wave-like nature of this 
observation (see figure [^, points to a mechanism with fixed pe¬ 
riodicity in time or space. This periodic structure, in combination 
with the detection of a strong velocity gradient resulting in ve- 
locity differences up to 7 x 10^ -9x10^ cm s“^ leads [Berne 


[et al.[ ([2010[) to propose that these ripples are manifestations of 
the KHI. 


Because of the high research interest in the Orion nebula and the 
surroundings regions, the physical conditions in the neighbour¬ 
hood of the observed ripples are fairly well documented, provid¬ 
ing an ideal case to numerically model the observed system. In 
[Berne & Matsumot^ ( [2012| ) an effort was undertaken to numeri¬ 
cally study the linear growth phase of a KHI with physical values 
deduced from observations. It was found that the used setup was 
indeed Kelvin-Helmholtz unstable for setups with magnetic field 
orientations close to perpendicular to the flow, and parallel to the 
separation layer between the Hjj and cloud region. 

In this work, our goal is to expand the numerical modelling of 
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Fig. 1. Observation of the ripples in Orion at 8 pm, taken with the 
Spitzer Infrared Array Camera. The spatial wavelength A, the orienta¬ 
tion of the phase velocity and the line ar regime length are iden - 
tified in the image. Credit: figure (1) from |B erne & Matsumoto| ( [2012| ), 
reproduced by permission of the AAS. 


the ripples in Orion in a way in which the observations can be 
directly compared to the modelling itself. To do so, several in¬ 
gredients are needed. First, the proposed setup (see sections |2.1| 
and |2.2| ) is simulated using a 3D numerical hydrodynamical sim¬ 
ulation from the start of the instability, through the linear phase 
and into the nonlinear phase. To perform these simulations we 
use th e MPI-AMRVAC code ( [Keppens et al.||2QT^ [Forth et aT 


2014), with numerical properties as described in section 2.3 In 


the mid-infrared observation a significant part of the radiation is 
due to dust emission. Therefore we use the gas-fdust module of 
the MPI-AMRVAC code to model the dynamics of dust particles, 
which are drag-coupled to the gas. We use a range of dust sizes 
and model it self-consistently with the gas dynamics. Finally, to 
connect the dynamical simulations to the observations we use 
the SKIRT dust radiative transfer code ( jBaes et al.|2Qll[[Camps| 
|& Baes|2015| ) to emulate the radiation by the dust particles and 
the effect of the ac tual geometry of the observed system, as ex¬ 
plained in section |2.4[ The properties of the outcome of these 
simulations are described in section [3] and the conclusions are 
discussed in section (4] 


2. Model 

2.1. Physical setup 

The setup used here is similar to that of the 2D setup of jBemej 
& Matsumotoj ( |2012| ), but here adjusted to a full 3D configura¬ 
tion. The domain of the simulation is a cube with L = 0.33 pc 
sides, and is initially divided in three regions along the y-axis: 
the upper part corresponds to the hot, low density H// region 
(njj = 3.34 X 10“^^ g cm"^, T// = 10"^ K), the lower part repre¬ 
sents the cold, high density molecular cloud (n^ = 1.67 x 10“^^ g 
cm“^, Tc = 20 K) and both are separated by a thin middle layer 
with thickness D = 0.01 pc. This boundary layer is thus oriented 
perpendicular to the y-axis. Note that the choice of density and 
temperature result in thermal pressure equilibrium between the 
upper and lower region as 


p=p 


hT 

mnp 


( 1 ) 


with p the pressure, the Boltzmann constant, Mh the mass of 
hydrogen and p the average molecular weight, set to yu = 1 here. 
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The energy density of the gas , e, can be calculated using the 
equation of state, and gives 


P Pv 
- 12 ’ 


with y = 5/3 the adiabatic constant and v the velocity of the flow. 


To initialise the dust content in the simulation domain, we 
assume that the d ust-to-gas mas s density ratio has the canoni¬ 
cal value of 0.01 ( |Spitzer||1954| ) in the molecular cloud region, 
and no dust is present in the hot H// region. We assume that 
the size distribution of dust particles, n, can be approximated 
as n{a) oc with the size of the particles, a, between 5 nm 
and 250 nm as was determined from excitation in the interstellar 
medium (ISM) by [Kim et al.| ( |1994) . We use four dust fluids to 
represent this power law size distribution with each fluid repre¬ 
senting a part of the size distribution, chosen in a way in which 
the total dust m ass in each dust fluid is the same (see [Hendrix &| 
[Keppens p014[ )). In this way, the resulting representative size of 
dust grain in the four dust fluids are 7.9 nm, 44.2 nm, 105 nm, 
and 189 nm, respectively. The grain densi ty of all dust fluids i s 
set to that of silicate grains, i.e. 3.3 g cm“^ (Draine & Lee 1984). 


The Hi I region has an initially uniform velocity of magnitude 


Vo 


10^ cm s Mn the direction parallel to our x-axis. 


Berne & 


[Matsumolo ( [2012[ ) propose that this high velocity is due to cham 


pagne flow, the resulting high velocity flow when the expanding 
H// breaks trough the molecular cloud. This velocity is similar 
to the shear velocity derived from observation in [Berne 
( [2010| ). In the molecular cloud region the velocity is initially set 
to zero. In contrast to [Berne & Matsumoto[ ( [2012| ), where a hy¬ 
perbolic tangent profile is used for both velocity and density, we 
use a linear profile in the middle layer that continuously links 
up with the constant velocities and densities on both si des of the 
layer. This is done in analogy with our previous work ( [Hendrix 


[& Keppen^[2014[ ), as it allows to better quantify the linear sta¬ 
bility properties. 

A perturbation is added by introducing an initial velocity com¬ 
ponent perpendicular to the boundary layer: 


/ ^ ,f>-3 I (y-My? \ . 

Vyp(x, y,z) =10 Vo exp--— [ sm (k^x) 


la] 


Icrl 


y 


-hlO %rect( —)(1 -2rand()), 


( 3 ) 


with (jy - 5D, cTz = L/5 and My and being the y- and z- 
coordinates of the middle point of the separation layer. The first 
part on the right side of equation ^ adds a sine perturbation with 
wavelength A = kxlln. We adopt T = 0.11 pc in accord with the 
observations in [Bgme et al.[ ( [2010| ). The second part on the right 
side of equation ^ adds random velocitie^ between -10“"^vo 
and 10“"^Vo in a layer of thickness 5D around the middle of the 
separation layer. The velocity in the z-direction is seeded with a 
similar random term: 


Vz,o(^, y,z) = 10 '‘vorectC^Xl - 2rand()). (4) 

The purpose of the exponential part in equation 0 in the y- 
direction is to preferentially locate the perturbation around the 
middle layer. The exponential part in the z-direction centres the 

^ The random function rand generates a random floating point value 
between 0 and 1, while the rect function (also called "rectangular func¬ 
tion") is one between -0.5 and 0.5 and zero elsewhere. 
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perturbation around the middle of the z-axis to confine the in¬ 
stability development region. These random perturbations in the 
velocity break the symmetry of the setup, and allow in essence 
all unstable modes to develop spontaneously, although the fixed 
A wavelength in the x-direction gets preference. 


2.2. Magnetic pressure 


Berne & Matsumot^ ( |2012| ) take into account a magnetic contri¬ 


bution in their 2D setup as well, assuming a uniform magnetic 
field with a strength of B = 200 pG in the entire domain based 
on observations of surrounding regions ( [Abel et JT]|2004[ [Bro¬ 
gan et al.|2005| ). Using the values of the physical setup (section 


2.1|) this results in a ratio between thermal and magnetic pressure 
Ppi = PtIPM = 0.0173, with Ppi the plasma beta value, meaning 
that the magnetic pressure is dominant over the thermal pressure 
contribution. The dominance of magnetic over thermal pressure 
is confirmed by observations in the orion molecular cloud (Berne 
et al.||2014|) , both for large and small scale structures. [Berne &' 


Matsumotq ( [20 1 2) note that the setup is most unstable when the 


magnetic field is perpendicular to the fiow and parallel to the 
contact layer. In this configuration, a uniform magnetic field only 
contributes as an additional magnetic pressure 




(5) 


This means that one can actually substitute the full MHD treat¬ 
ment by a HD treatment with an additional pressure term, in 
which the total pressure is raised while keeping the density fixed 
(thus artificially increasing the temperature). When calculating 
the thermal energy of the gas to quantify the coupling to the dust 
(see ( [Forth et al.[2014[ )), this artificial term is subtracted to obtain 
the relevant temperature. To demonstrate that this approximation 
is valid, we compare evolution of an MHD setup with that of a 
HD -r Pm simulation in section [TT] 


2.3. Numerical method 


We us e the MPI-AMRVAC code ( [Keppens et al.|2012} [Forth et~aL| 
2014[ ) for all the hydrodynamical (HD) and magnetohydrody- 
namical (MHD) sim ulations. The dust module of MPI-AMRVAC, 
discussed in detail in [Hendrix & Keppens [ ( [20 14| ), allows to add 
dust to a HD simulation by adding multiple dust fiuids. These 
fluids follow the Euler equations with vanishing pressure (LeV- 
eque[2004] ) and couple to the gas fluid through a drag force term 
Each dust fluid has its own physical properties such as grain size 
and grain material density. Typically we use multiple dust fluids 
with the same grain material density and different grain sizes to 
model the size distribution in the ISM. 

Eor the 3D simulations we use four levels of adaptive mesh re¬ 
finement (AMR), resulting in an effective resolution of 448 x 
1792 X 448 cells. The triggering of extra refinement levels is 
based on a combination of the gradients in the gas fluid and 
those in the dust fluid representing the largest grains. Because 
the actual physical domain is cube shaped, this resolution results 
in a four time higher resolution perpendicular to the flow (see 
section mi- This is necessary to resolve all small-scale varia¬ 
tions that develop during the linear (and also the nonlinear) phase 
of the instability. The solution of the coupled gas-fdust fluid 
equations is advanced using a total variation diminishing Lax- 
Eriedrich (TVDLE) scheme with a two-step predictor-corrector 
time discretisation and a monotonised central (MC) type limiter 
( van Leer|1977[ ). To ensure stable time-stepping the timestep is 



Fig. 2. Growth of the kinetic energy perpendicular to the bulk flow. The 
MHD and HD simulation that take into account the magnetic pressure 
are similar, while the HD simulation without magnetic pressure behaves 
differently. The 3D setup is also shown up to t = 0.01 and has a growth 
rate similar to that of the 2D setup. 


limited by using a CEL number of 0.6 for gas and dust, as well 
a separate dust acceleration criterion based on the stopping time 
of dust grains ( [Laibe & Frice|2012] ). 

2.4. Radiative transfer 

To be able to directly compare the output from the 3D hydro- 
dynamical simulations with observations, post-processing of the 
data is performed with the Monte Carlo radiative transfer code 
SKIRT ( [Baes et al.[|WT| [Camps & Bae^[MT5] ). SKIRT sim¬ 
ulates continuum radiation transfer in dusty astrophysical sys¬ 
tems by launching a set of photon packages in a given wave¬ 
length range through the dust distribution obtained from our 
dynamical simulations. These packages are followed for sev¬ 
eral cycles of multiple anisotropic scattering, absorption and (re- 
)emission by interstellar dust, including non-local thermal equi¬ 
librium dust emission by transiently heated small grains. Emis¬ 
sion from stochastically heated grains is used in all the results in 
this work and typically around 4 dust emission cycles are needed 
to come to equilibrium. 

To launch the packages into the domain, we use a (stellar) point- 
source at a given distance outside of the simulated domain as 
our source of initial photons. Fhoton packages in a wavelength 
range between 0.01 pm and 1000 pm are incorporated. In SKIRT 
we use exactly the same distribution of dust species as the one 
obtained from MFI-AMRVAC, meaning that the mass density 
distribution of the four dust fluids is used for each representa¬ 
tive part of the grain size distribution and that, just like in the 
HD simulations, we adopt silicate properties for the grains in the 
radiative transfer. 


3. Results 

3.1. 2D analysis 

To prove that an MHD setup with the magnetic field component 
perpendicular to the flow direction and parallel to the bound¬ 
ary layer can be reasonably approximated by a similar setup in 
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Fig. 3. Gas density plots of the KHI in 2D and 2.5D after the end of the 
linear phase. The density units are in g cm"^. In all figures the entire do¬ 
main (0.33 pc X 0.33 pc) is shown. Left: A 2D simulation of the KHI in 
HD with dust and an artificial magnetic pressure term pm added to the 
total pressure at ^ = 0.007 (6.84x10"^ years). Centre: The same setup, 
but in 2.5D MHD with a magnetic field perpendicular to the plane, also 
at ? = 0.007. Right: A 2D HD simulation without the effect of a mag¬ 
netic field added into the total gas pressure, at ? = 0.02 (1.95x10^ years). 
Note this figure is taken at a different time as the linear phase end later 
in this case. 


HD but with add ed pr essure, we simulate the setup discussed in 
sections [TT] and first in 2D, but in three variations: a HD 
simulation without a magnetic contribution, an MHD setup with 
magnetic field, and an HD simulation with the magnetic field 
contribution added to the pressure. The MHD setup is actually 
simulated in 2.5D, as it includes the information of the veloc¬ 
ity and magnetic field perpendicular to the simulated plane. The 
simulated plane in 2D corresponds to a slice in the 3D simula¬ 
tion perpendicular to the v-y plane and through the centre of the 
simulated domain. In figure [2 the buildup of kinetic energy per¬ 
pendicular to the flow direction is shown for all three 2D setups, 
and for the 3D run discussed further on. Clearly, for the MHD 
setup and the HD plus magnetic pressure setup the growth rate 
in the linear regime (up iot = 0.006 in code units, or ~ 5.87x10'^ 
years) is the same. The growth rate is significantly slower when 
the magnetic pressure is ignored. Also, figure shows that the 
formed structures are of similar size and shape in the two simu¬ 
lations where the magnetic pressure is taken into account. Small 
differences include the formation of small-scale structures on top 
of the larger structure. These small-scale perturbations are also 
present in the HD setup, but develop faster in the MHD simu¬ 
lation. The reason that they are less apparent in the HD simula¬ 
tion is because in the MHD case they seemingly grow faster due 
to small inhomogeneities (a decrease by 2%) in the magnetic 
field, leading to numerical differences that accumulate over time. 
When the magnetic pressure is not taken into account, it can be 
seen in figure that the morphology is very different. Because 
the total pressure is lower, the Mach number for the flow at the 
boundary is higher, causing shocks to propagate. These shocks 
also cause the striped structure in the high density region. We 
will now further discuss a full 3D gas plus dust setup that has the 
pressure adjusted to account for the magnetic pressure effects. 


3.2. 3D model 

In figurej^it can be seen that the growth rate of the 3D simulation 
is comparable to that of the 2D simulations in which the effect 
of the magnetic field is taken into account. Due to the added 
computational cost in 3D, this simulation is only followed until 
r = 0.01 in code units, or up to about 9.78x10"^ year. 
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time 

Fig. 4. Time evolution of the maximal density enhancements in the 3D 
simulation for all four dust fluids, with dust 1 representing the smallest 
grains (7.9 nm) and dust 4 the largest grains (189 nm). 


3.2.1. Dust distribution 


In previous work ( Hendrix & Keppens|2014| ) we found that in a 
3D setup with the same density on both sides of the separation 
layer, the KHI can cause the dust density to increase by almost 
two orders of magnitude. These strong increases in dust density 
occur in filament-like locations between the vortices when 
dust is swirled out of the vortices and compressed into these 
regions. This process if strengthened further by additional 3D 
instabilities. Also, it was found that the process of dust density 
enhancement is stronger for larger dust particle sizes. Figure 
shows that in the setup used here the growth in local dust 
density is less strong. During the end of the linear phase, i.e. 
up to time t = 0.006 in figure]^ the maximal density increases 
gradually, and the rate of increase is proportional to the grain 
size. In the further nonlinear stage the densities still increase, 
however the relation between instantaneous local maximal 
density and grain size gets modified. Similarly to what was seen 
in [Hendrix & Keppens| ( [?014 ), the density enhancements are 
significantly stronger in 3D than in 2D, where the maximum 
increase is less than 15% for all dust species in the 2D case with 
magnetic pressure added. Clearly, 3D effects are paramount 
when studying dust growth. 


The dust density enhancements are strongest in three distinct re¬ 
gions, which are indicated in figure Chronologically dust first 
accumulates in the convex outer region of the KH wave (the re¬ 
gion labeled with 1 in figure [^. This is due to the acceleration 
of dust by gas in the concave region when the gas swirls around 
the low pressure region created by the KHI. Next, the arc-like 
structure below the surface of the wave, i.e. region number 2 in 
figure is formed. This region forms when the KHI accelerates 
the bulk of the gas upward into the low density region, and the 
dust is dragged with it. The location of the region is caused by a 
gradient in the drag strength, as the velocity difference between 
gas and dust is stronger under the region than above, causing the 
underlying dust to overtake the dust above it. The third dust gath¬ 
ering region is along the boundary between high and low density 
regions in between two successive waves or KHI rolls. A dust 
pile-up is seen here in the nonlinear stage when the velocity of 
the gas around the low pressure vortex is highest. In animated 
views one can see how the end point of the flow that passes over 
the crest of the waves moves from location 1 to a spread out 
region all along the density boundary, i.e. up to location 3 as in¬ 
dicated. 
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Fig. 5. Density of the largest dust species {a - 189 nm) in a slice from 
the 3D simulation (z = 0.165pc) at t = 0.0065 years). Only a 

part of the simulated region with an extend of 0.138 pc in the v-direction 
is shown. Three distinct regions of dust density enhancement are indi¬ 
cated with labels 1, 2 and 3 discussed in the text. The velocity field of 
the largest dust species in the x - y plane is indicted with the use of 
vectors, the largest velocity are around 6 x 10^ cm s"^ 



Fig. 6. Volume plot of the total dust density at t = 0.01 (9.78x10"^ years). 
Only densities higher than the initial maximum density {pd - 1.67 X 
10"^^ g cm"^) are visualised. 


While dust density increases up to a factor 10 are observed in 
these three regions for the four dust species, the actual location 
of these dust-gathering regions does not necessarily fully coin¬ 


cide for all dust species, similar to the findings in [Hendrix & 
|Keppens| ( [2014| ) where a clear size-separation was evident. Also, 
the actual importance of the three regions is distinct for different 
grain sizes. Therefore, the increase of the total dust density will 
be less strong and distributed over a larger region. Furthermore, 
the strongest increases can be found in small local clumps, as 
can be seen in figure visualising the total dust density con¬ 
centrations. Quantitatively speaking, while 14.76% of the total 
volume experiences a total dust density enhancement of more 
than 5%, in only 0.03% of the total volume the total dust den¬ 
sity more than doubles (regions indicated in orange a nd red in 
figure [6|). Th is is in contrast with the 3D simulations in |Hendri^ 
|& Keppen^ ( |2014| ), where the high density dust is found in long 
filamentary structures and more than 4.5 % of the volume ex¬ 
hibits a doubling of the total dust density. The main differences 
reside in the adopted initial density contrast, as well as the fact 
that here only the molecular cloud region initially had dust. 



Fig. 7. Geometry of the stellar object (photon source) and observer loca¬ 
tion with respect to the structures in Orion, designated by independent 
angles a and p, respectively. In this image, the location of the source 
and observer are shown with respect to the KH features at t=0.084 
(8.21 xlO"^ years). The black-white image is actually a SKIRT image 
at 54 /rm, where we see the radiation which is coming from dense and 
heated dust in the billow structures formed by the KHI. In this image, 
the observer is located perpendicular to the x - y plane. 


3.3. Modelling observations 

In the previous section we have outlined how the model setup 
from section |2.1| evolves into a nonlinear 3D KHI. Next, we 
investigate how the simulated structures would look in synthetic 
observations. As described in section [2^ the dust distribution 
of our 3D simulations is used as input for the SKIRT radiative 
transfer code. To see to which degree our simulations corre¬ 
spond to the actual observed structures (figure [^, in addition 
to the hydrodynamical setup one has to take into account the 
orientation in relat ion to the observer, as well as the location of 
the light source(s). Berne et ^ ( [2010 ) indicated that the star 6^ 
Orionis C, a massive type 07V star ( Donati et al.||2002{ |Wade 
et al.||2006| ) located in the H// Trapezium region at a distance 


of ~ 3.4 pc from the cloud, illuminates the ripples from behind 
with respect to the observer. In SKIRT the radiation of this star 
is simulated by adding a point source of photons at J = 3.4 pc 
and inclination a with respect to the initial separation layer in 
the HD simulation, as illustrated in figur e [7| For the ra diation 
of the star we use a model spectrum from [Martins et al.|p005| ) 
with corresponds to a star with physical properties comparable 
to those of 6^ Orionis ca The location of the observer with 
respect to the simulated domain must also be specified in 
SKIRT. As shown in figure |7] the observer is placed at an angle 
P with respect to the initial separation layer in the HD simulation. 


Because the actual inclination between the observer, the 
billows and the background radiation source are hard to gauge 
from the observation, several different values of a and p were 
tried to investigate their role. Table gives an overview of 
several SKIRT geometries we will discuss here. An interesting 
setup to look at first is case D (figure top right). With this 
arbitrary choice for the geometry {a = 60° and p = 90°) the 
result is rather different from the observations. While some 
periodicity is observable, no sharp elongated structures are seen. 
The diffuseness of the radiation in case D can be seen to be 

^ Model T46pl_logg4p05.sed from http://www.mpe.mpg.de/ 
~martins/SED.html 
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Table 1. Summary of the SKIRT radiative transfer models, with a the 
angle between the star and the cloud, and (5 the angle between the cloud 
and the observer (see figurej^, and the time in code units. 


Case 

a 

p 

time 

A 

40 

128 

0.0082 

B 

40 

90 

0.0082 

C 

60 

128 

0.0082 

D 

60 

90 

0.0082 

E 

51 

90 

0.0082 

F 

51 

128 

0.0082 

G 

51 

90 

0.01 

H 

51 

128 

0.01 


inherent to an observer angle of 90°. Figure demonstrates 
that when going from t = 0.0082 in E to r = 0.01 in G, 
while the onset of the nonlinear phase increases the d evelop - 
ment of small-scale features (as discussed in section |3.2.1| ), 
the emission in the nonlinear phase remains diffuse in both cases. 

In figure [7] we see that the emission at 54 pm is strongest 
where the dust is directly radiated by the source, but the colder 
dust inside the KH billows also radiates at this wavelength. At 
shorter wavelengths such as 8.25 pm, the direct light is the more 
important and only dust close to the edges of the billows radi¬ 
ates. To get features more reminiscent of the observations we 
can use this knowledge to consider two changes to the geome¬ 
try of the source and the observer. On the one hand, the angle a 
can be chosen to maximise the photons from the source reach¬ 
ing the protruding billows and not the rest of the cloud, which 
increases the amount of observed photons in a more compact lo¬ 
cation. Nevertheless, the effect of changing a is small at 8.25 pm, 
as demonstrated by comparing cases A to C and B to D in figure 
On the other hand the observers angle p can be chosen to be 
along the billows, maximising the perceived compactness. The 
change in observer angle has a much stronger impact. Changing 
P from 90° in case B to yS = 128° in case A clearly decreases 
the thickness of the features, increases the flux in the elongated 
regions, and enhances the contrast between the bright en dark 
regions. The choice for “optimal angles" is illustrated in figure 
The values we find are a = 5\° and p = 128°. These val¬ 
ues are used in cases F and H (figure [T^. Using this geometry, 
a fair approximation of the real observations can be made, at a 
comparable wavelength. The evolution from case F into H again 
displays the formation of the small scale structures in the non¬ 
linear phase, on a scale which is comparable to the local bends 
in the infrared observations. 


4. Conclusions 

In the previous sections, we have modelled a region of the Orion 
molecular cloud in which elongated ripple features are observed. 
To do so, we have built upon previous numerical models, and 
expanded these to full 3D dusty hydrodynamics coupled to a 
radiation transfer code designed for simulating dusty astrophys- 
ical systems. The synthetic images allow a direct comparison 
with the observations. In the infrared observations, the ripples 
are thin, elongated features that have a clear periodicity and 
are sharp and bright compared to the background radiation. All 
these features can also be reproduced by our model. The hy- 
drodynamical simulations confirm that the previously proposed 
setup is indeed KH unstable for the observed spatial wavelength. 
We find that the dynamical contribution of dust with a size 
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Fig. 8. SKIRT simulations of the same dataset with different geome¬ 
tries. From left to right and top to bottom: B, D, A, C. Horizontally the 
observers angle p is the same {p - 90° on iop,p - 128° below) and the 
same scaling is used. Note that the flux quantification is arbitrary here 
and no effort has been taken to compare these to real values. Vertically 
the irradiation angle is constant (a = 40° left, a = 60° right). All images 
are observed at 8.25 pm. 



Fig. 9. Synthetic observation of the KHI at 8.25 pm, with fixed observa¬ 
tional angle p - 90° and a - 128° (cases E and G). Two different times 
are shown, left: t - 0.0084, right: t - 0.01 or 8.21x10"^ and 9.78x10"^ 
year, respectively). During this interval the development of small-scale 
perturbations in the nonlinear phase can be seen. A linear scale is used 
for the intensity of the images. 


distribution typical for the ISM does not inhibit the formation of 
the KHI, and the growth rate in 3D is similar to that of the 2D 
simulation. We see that the presence of a background star is able 
to light up the features of the KH billows. Also, the synthetic 
images demonstrate clearly that the geometry is of great impor¬ 
tance in distinguishing the KH features from the background. 
Observers located in a direction perpendicular to the shearing 
layer would observe some periodicity, however with shallow 
features over a continuous background, while observers which 
look along the formed billows observe them very sharp and 
bright compared to the background. Nevertheless, even when 
considering the most optimal geometry, the ripples are still 
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Fig. 10. Synthetic observation of the KHI at 8.25 yum, with observational 
angle p - 128° and cr = 51° (cases F and H). Two different times are 
shown, left: t - 0.0084, right: t - 0.01 or 8.21x10"^ and 9.78x10"^ 
year, respectively). In comparison to the images at/5 = 90°, the features 
of the KHI are more pronounced and clearly distinguishable from the 
background. A linear scale is used for the intensity of the images. 


somewhat wider than the sharp ripples of the observations. 
Additional to geometrical effects, the sharp features may point 
to strong local density increases in the dust , however in contrast 
to our previous investigation of dusty KHI ( [Hendrix & Keppeii^ 
|2QI4| ) only small increases in dust density are seen here, and the 
highest increases are found in small and compact clumps and 
not elongated regions. The treatment of additional physics such 
as self gravity and magnetic fields may lead to these additional 
density increases as was shown for larger scale structures in jVanj 
[Loo et al!] ( |2014| ). It is unclear if a significant effect would also 
be expected here, as in sectionjSHjthe magnetic field only causes 
minor deviations in the 2D setup. For simulations in 3D, the 
strong magnetic field (plasma= 0.0173) may somewhat alter 
the outcome of the simulations in the nonlinear phase, when 
secondary 3D instabilities break the earlier quasi-2D behaviour. 
|Ryu et al.j ( |2000| ) demonstrated that even weak magnetic fields 
can be of importance in the nonlinear regime. While a strong 
magnetic fields may suppress the growth of hydrodynamical 
perturbations perpendicular to the fields, [Matsumoto & Sela| 
(2007]) find that in cases with plasma beta as low as Ppi = 0.1 
secondary 3D instabilities also occur and cause small scale 
fragmentation along the initial magnetic field, however at a stage 
far in the nonlinear regime. The resulting influence of the 3D 
magnetic field on the dynamics of the dust grains, and thus also 
the observed structures, is further complicated by the unknown 
charge of the dust grains. While for example [Hoang et aT] ( |2012| ) 
have calculated mean grain charging as function of grain sizes 
for different ISM phases, the charging of grains can be location 
dependant due to for example interaction with a radiation field, 
as is the case here. Fully taking into account the magnetic 
field would thus also require further assumptions to be made 
with regard to dust distribution as a function of the both the 
size and the charge. Furthermore, the strength of the magnetic 
field is one of the less constrained parameters in the model; 
while the value in the model = 200 yuG) is representative for 
surrounding regions, no local measurements of orientation and 
strength exist to our knowledge. As the magnetic pressure is 
shown to be of importance in finding the correct value for the 
growth rate (section 113, the outcome would be different if a 
different magnetic field was assumed. This would especially be 
the case for different relative orientations of this field and the 
flow shear. 


Another important factor which may change the outcome 
of the simulations is the actual width of the shearing layer 


between the hot medium and the molecular cloud. The width 
is an important parameter in the evaluation of the stability and 
growth of the KHI instability. The 
pc) is in analogy with the value of 
where it is argued that this value represents the width of 
the photodissociation region (PDR), where molecular gas is 
dissociated by the far ultraviolet photons of the background star 
0^ Orionis C. Nevertheless, as discussed in the supplement of 
Berne et J^p0I0[ ), actually a broader (~ 0.1 pc) photo-ablation 
region forms between the PDR and the hot medium. Due to 
its thickness this region may inhibit the formation of the KHI 
with wavelengths in the range of the observed periodicity in the 
ripples or shorter, as a boundary layer of thickness D inhibits the 
growth of perturbations with A < 4.9ID ( [Chandrasekhar! 1961 [ 
[Hendrix & Keppens|2014[ ). Additionally it should be noted that 
the effect of heat conduction, which has not been included in 
this work, can be of importance in the formation of the shearing 
layer between the hot me dium and the molecular cloud. Indeed, 
[Vieser & Hensler[ ( [2007 [ ) demonstrate that heat conduction can 
reduce the steepness of the velocity gradient between the cloud 
and a streaming flow, stabilising the surface of the cloud against 
the development of the KHI. 

While these remarks demonstrate that additional physics 
may be needed to understand the full range of interactions oc¬ 
curring in the Orion nebula, in this work we tried to model the 
observations of its KH ripples in full detail. We demonstrated 
that a full treatment of gas and dust dynamics, including a range 
of dust sizes, coupled with radiative transfer provides a promis¬ 
ing approach to explaining the observations. Even though the 
physical values in the models are prone to intrinsic observational 
uncertainties or assumptions, we see that these values are reason¬ 
able in reproducing most of the features when the most optimal 
geometrical model is used. 
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